---
title: "replication_testingH2"
author: "Myung Jung Kim"
date: "2026-01-30"
output: html_document
---


```{r packages load}
library(readr)
library(tidyverse)
library(marginaleffects)
library(sandwich)
library(graphics)
library(stats)
library(ggplot2)
library(ggeffects)
library(sandwich)
library(lmtest)
```


```{r H2: Analysis to test H2 (ICC) & Table 5}
library(readr)
df <- read_csv("df2_triadyear_amnesty.csv")

names(df)

df <- df%>%dplyr::select(dyad_old_id,dyad_new_id, side_a, side_a_id, side_a_id_old,
                   side_b, side_b_id, side_b_id_old, year, conflictid, country, ccode,time,
                   dyad_yearatwar, dyad_yearatwarln, numdyads, intensitylevel, TBR_NSAdummy, 
                   rebcap,blood,
                   start_year,end_year,start_post98,start_post2002,end_pre98,end_pre2002,
                   v2x_rule_Home,Home_v2svstterr,
                   hram, 
                   Host_mean_v2svstterr, 
                   rebel_OSV_dummy, deaths_civilians, osv_reb_cumulative,
                   PKO_dummy_Home,
                   hram, dhram, dhram1, thram, tnhram, presname, 
                   Host1,Host2,Host3,Host4,Host5,Host6,
                   Host1_COW, Host2_COW,Host3_COW, Host4_COW,Host5_COW,Host6_COW,
                   ICC_Involve_lniccactionscum_cl2_Home, 
                   ICC_Involve_lniccactionscum_cl2_Host1,
                   ICC_Involve_lniccactionscum_cl2_Host2,
                   ICC_Involve_lniccactionscum_cl2_Host3,
                   ICC_Involve_lniccactionscum_cl2_Host4,
                   ICC_Involve_lniccactionscum_cl2_Host5,
                   ICC_Involve_lniccactionscum_cl2_Host6,
                   Host_all_mean_ICC_Involve_lniccactionscum_cl2,
                   PKO_dummy_Host1, PKO_dummy_Host2,PKO_dummy_Host3,PKO_dummy_Host4,
                   PKO_dummy_Host5,PKO_dummy_Host6,PKO_Host_all_dum,
                   Host1_v2svstterr,Host2_v2svstterr,Host3_v2svstterr,
                   Host4_v2svstterr, Host5_v2svstterr,Host6_v2svstterr, 
                   Host_mean_v2svstterr, post98,post2002)

dim(df)
# cannot find original df_NSA so trying this: 
df_NSA <- df%>%filter(TBR_NSAdummy ==1)
df_NSA_post2002<- df_NSA %>%filter(year>2001) #from when ICC involvement was possible

table(df_NSA_post2002$Host_all_mean_ICC_Involve_lniccactionscum_cl2, df_NSA_post2002$year)

df_NSA_post2002$dyad_yearatwarln <- log1p(df_NSA_post2002$time)


## ################## Host ICC (mean)   #########################
AX_icc_gradual1 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2  
         + v2x_rule_Home 
         + Home_v2svstterr
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual1_cluster <- coeftest(AX_icc_gradual1, 
                                    vcov = sandwich::vcovCL(AX_icc_gradual1, 
                                                            cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual1_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual2 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2  
         + v2x_rule_Home
         + Home_v2svstterr
         +Host_mean_v2svstterr
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual2_cluster <- coeftest(AX_icc_gradual2, vcov = sandwich::vcovCL(AX_icc_gradual2, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual2_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual3 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2  
        + v2x_rule_Home
         + Home_v2svstterr
         +Host_mean_v2svstterr
         + numdyads 
        +intensitylevel 
         +dyad_yearatwarln 
         ,data=df_NSA_post2002, family = binomial)


AX_icc_gradual3_cluster <- coeftest(AX_icc_gradual3, vcov = sandwich::vcovCL(AX_icc_gradual3, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual3_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual4 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2  
     + v2x_rule_Home
         + Home_v2svstterr
         +Host_mean_v2svstterr
         + numdyads 
        +intensitylevel 
         +dyad_yearatwarln 
         + rebel_OSV_dummy
         +rebcap
         + blood
     +PKO_dummy_Home
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual4_cluster <- coeftest(AX_icc_gradual4, vcov = sandwich::vcovCL(AX_icc_gradual4, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual4_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual5 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2 
     + v2x_rule_Home
         + Home_v2svstterr
         +Host_mean_v2svstterr
         + numdyads 
        +intensitylevel 
         +dyad_yearatwarln 
          + rebel_OSV_dummy
              +rebcap
         + blood
     +PKO_dummy_Home
     +PKO_Host_all_dum
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual5_cluster <- coeftest(AX_icc_gradual5, vcov = sandwich::vcovCL(AX_icc_gradual5, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual5_cluster


library(stargazer)
library(texreg)
# Using texreg to display the models
texreg(list(AX_icc_gradual1_cluster,AX_icc_gradual2_cluster,AX_icc_gradual3_cluster,AX_icc_gradual4_cluster,AX_icc_gradual5_cluster), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

# For numbers and AIC info, See this: 
texreg(list(AX_icc_gradual1,AX_icc_gradual2,AX_icc_gradual3,AX_icc_gradual4,AX_icc_gradual5), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

## figure

library(ggeffects)
library(ggplot2)
library(sandwich)

# predicted probability for ICC Involvement (discussed in paper)
vcov_matrix <- vcovCL(AX_icc_gradual5, cluster = df_NSA_post2002$side_a_id)

pred_icc_final <- ggpredict(AX_icc_gradual5, 
                            terms = "Host_all_mean_ICC_Involve_lniccactionscum_cl2 [all]",
                            vcov.fun = vcov_matrix)

pred_icc_final
```

```{r H2: Appendix - Cumulative ICC analysis, Table 5}
## ################## Host ICC (mean)   #########################
AX_icc_gradual1 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2
        + v2x_rule_Home 
         + Home_v2svstterr
         +intensitylevel 
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual1_cluster <- coeftest(AX_icc_gradual1, vcov = sandwich::vcovCL(AX_icc_gradual1, cluster =df_NSA_post2002$side_a, type = "HC0"))

AX_icc_gradual1_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual2 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2 
         + v2x_rule_Home
         + Home_v2svstterr
         +intensitylevel 
         +Host_mean_v2svstterr
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual2_cluster <- coeftest(AX_icc_gradual2, vcov = sandwich::vcovCL(AX_icc_gradual2, cluster =df_NSA_post2002$side_a, type = "HC0"))

AX_icc_gradual2_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual3 <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2
        + v2x_rule_Home
         + Home_v2svstterr
         +intensitylevel 
         +Host_mean_v2svstterr
        + numdyads 
        +dyad_yearatwarln 
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual3_cluster <- coeftest(AX_icc_gradual3, vcov = sandwich::vcovCL(AX_icc_gradual3, cluster =df_NSA_post2002$side_a, type = "HC0"))

AX_icc_gradual3_cluster

## ################## Host ICC (mean)   #########################
AX_icc_gradual4_cumosv <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2  
      + v2x_rule_Home
         + Home_v2svstterr
         +intensitylevel 
         +Host_mean_v2svstterr
         + numdyads 
         +dyad_yearatwarln 
      +osv_reb_cumulative 
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual4_cluster_cumosv <- coeftest(AX_icc_gradual4_cumosv, vcov = sandwich::vcovCL(AX_icc_gradual4_cumosv, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual4_cluster


## ################## Host ICC (mean)   #########################
AX_icc_gradual5_cumosv <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2
     + v2x_rule_Home
         + Home_v2svstterr
         +Host_mean_v2svstterr
        +dyad_yearatwarln 
         +intensitylevel 
        + numdyads 
      +osv_reb_cumulative  
  +rebcap+ blood 
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual5_cluster_cumosv <- coeftest(AX_icc_gradual5_cumosv, vcov = sandwich::vcovCL(AX_icc_gradual5_cumosv, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual5_cluster_cumosv


## ################## Host ICC (mean)   #########################
AX_icc_gradual6_cumosv <- glm(hram ~ Host_all_mean_ICC_Involve_lniccactionscum_cl2 
   + v2x_rule_Home
         + Home_v2svstterr
         +Host_mean_v2svstterr
        +dyad_yearatwarln 
         +intensitylevel 
        + numdyads 
      +osv_reb_cumulative  
  +rebcap+ blood 
    +PKO_dummy_Home
      +PKO_Host_all_dum
         ,data=df_NSA_post2002, family = binomial)

AX_icc_gradual6_cluster_cumosv <- coeftest(AX_icc_gradual6_cumosv, vcov = sandwich::vcovCL(AX_icc_gradual6_cumosv, cluster =df_NSA_post2002$side_a_id, type = "HC0"))

AX_icc_gradual6_cluster_cumosv


library(stargazer)
library(texreg)
# Using texreg to display the models
texreg(list(AX_icc_gradual1_cluster,AX_icc_gradual2_cluster,AX_icc_gradual3_cluster,AX_icc_gradual4_cluster_cumosv,AX_icc_gradual5_cluster_cumosv,AX_icc_gradual6_cluster_cumosv), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

# For N and AIC info, see this: 
texreg(list(AX_icc_gradual1,AX_icc_gradual2,AX_icc_gradual3,AX_icc_gradual4_cumosv,AX_icc_gradual5_cumosv,AX_icc_gradual6_cumosv), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))
```

